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Triangulating  Trimmed  NURBS  Surfaces 

Chang  Shu  and  Pierre  Boulanger 


Abstract.  This  paper  describes  techniques  for  the  piecewise  linear  ap- 
proximation of  trimmed  NURBS  surfaces.  The  problem,  called  surface 
triangulation,  arises  from  many  applications  in  CAD  and  graphics.  The 
new  method  generates  triangular  meshes  that  are  adaptive  to  the  local 
surface  curvature.  We  use  efficient  data  structures  for  the  handling  of 
trimming  curves.  We  also  generate  Delaunay  triangulation  on  the  surface 
to  improve  the  quality  of  the  meshes. 


§1.  Introduction 

Tensor-product  NURBS  are  widely  used  in  today’s  CAD  systems  for  describing 
and  exchanging  surface  geometry.  For  many  applications,  however,  piecewise 
linear  approximations  of  smooth  surfaces  are  required.  Examples  of  these  ap- 
plications include  finite  element  analysis,  stereo-lithography,  and  visualization 
of  geometric  models.  In  these  applications,  we  need  to  generate  a triangular 
mesh  that  approximates  the  original  surface  within  a given  tolerance.  We  refer 
to  this  problem  as  surface  triangulation,  stated  in  the  following  definition. 

Definition.  Given  a NURBS  surface  N(u,  v),  its  trimming  boundary,  and 
a real  number  e,  the  surface  triangulation  problem  is  to  find  a set  of  linearly 
parameterized  triangles  {Tj}  such  that 

1)  Any  triangle  T*  satisfies  sup  ||Tj(w,  v)  — N(u,  v)||  < e. 

2)  For  any  triangle  edge  not  on  the  boundary,  there  is  exactly  one  neighbor- 
ing triangle  sharing  this  edge. 

The  first  condition  is  usually  called  chord  height  tolerance,  which  restricts 
triangles  to  be  close  to  the  surface.  The  second  condition  requires  the  trian- 
gular mesh  to  be  topologically  correct. 

A good  surface  triangulation  algorithm  is  expected  to  be  efficient  because 
real  world  models  tend  to  contain  large  numbers  of  surface  patches.  Further- 
more, certain  optimization  factors  are  desirable.  Two  of  the  most  important 
ones  are  triangle  shape  and  the  number  of  triangles  in  the  mesh.  For  exam- 
ple, in  finite  element  analysis,  triangles  with  bad  aspect  ratio  (one  angle  is 
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significantly  smaller  or  larger  than  the  others)  reduce  the  solution  precision. 
In  all  applications,  a mesh  with  a small  number  of  triangles  saves  computing 
and  transmission  time  as  well  as  storage  space. 

Several  authors  [4,6,8,14]  have  approached  the  surface  triangulation  prob- 
lem by  computing  a bound  on  the  length  of  triangle  edges  in  parametric  space 
so  that  if  all  triangles  have  their  edge  lengths  smaller  than  the  bound,  the 
resulting  triangulation  satisfies  the  chord  height  tolerance.  Since  the  edge 
length  bound  applies  to  the  whole  surface,  the  density  of  the  triangulation 
distributes  uniformly  across  the  surface  and  may  lead  to  unnecessarily  large 
mesh  size.  Along  another  line  of  thought,  Klein  and  Strafier  [5]  considered 
the  problem  of  placing  points  based  on  the  surface  curvature.  Recently,  Piegl 
and  Tiller  [11]  used  adaptive  subdivision  of  the  surface.  Obviously,  for  a given 
chord  height  tolerance,  adaptive  algorithms  generate  fewer  triangles  than  the 
uniform  subdivision  algorithms.  But  adaptive  algorithms  tend  to  be  slower. 

In  this  paper,  we  give  a method  that  has  the  following  features: 

1)  adaptive  to  the  surface  curvature, 

2)  efficient  insertion  of  the  trimming  curves, 

3)  triangle  shape  improvement. 

Our  general  strategy  is  that  we  first  approximate  the  surface  with  hierar- 
chical quadrilaterals  without  considering  the  trimming  curves,  then  we  insert 
the  trimming  curves  and  triangulate  the  quadrilaterals.  The  result  is  a trian- 
gulation that  satisfies  the  chord  height  tolerance.  We  improve  the  efficiency 
of  trimming  curve  insertion  by  organizing  the  quadrilateral  hierarchies  in  a 
quadtree  structure.  Also,  we  improve  the  quality  of  the  triangles  by  convert- 
ing the  initial  triangulation  to  a Delaunay  triangulation. 

§2.  Curve  and  Surface  Subdivision 

We  begin  by  discretizing  the  surface  and  its  trim  curves  independently.  We 
assume  that  the  surfaces  are  trimmed  in  the  parametric  domain  and  the  trim- 
ming curves  are  represented  as  NURBS  with  consistent  orientation.  Our  first 
objective  is  to  approximate  the  curves  with  connected  line  segments  such  that 
they  do  not  deviate  from  the  surface  more  than  the  tolerance  e.  From  well- 
known  results  in  B-Spline  theory  [7,10],  a NURBS  curve  can  be  split  into  two 
pieces  without  changing  its  shape  by  inserting  new  knots.  The  consequence 
of  this  splitting  is  that  we  introduce  new  control  points  that  are  closer  to  the 
curve  than  the  control  points  of  the  original  curve.  If  we  keep  dividing  in  this 
way,  the  control  polygon  converges  to  the  surface.  When  a sub-curve’s  control 
polygon  becomes  “flat”  enough,  we  can  stop  the  dividing  process.  Accord- 
ing to  the  convex  hull  property  of  NURBS,  the  maximum  distance  from  any 
point  on  the  sub-curve  to  the  line  segment  joining  the  two  end  control  points 
is  bounded  by  the  maximum  distance  between  control  points  to  the  segment. 
Therefore,  we  can  use  this  bound  to  control  the  flatness  of  the  sub-curves. 

The  surface  can  be  approximated  in  the  same  way  by  quadrilaterals.  At 
this  time,  we  ignore  the  trimming  boundary.  Here,  we  insert  knots  in  both  u 
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Fig.  1.  Full  surface  subdivision. 

and  v directions.  The  flatness  test  is  a little  more  complex.  We  examine  every 
row  and  column  of  the  control  polygon  and  test  their  flatness.  We  also  have 
to  consider  the  twist  factor  of  a surface  patch.  Peterson  [9]  gives  a subdivision 
method,  which  we  generally  follow.  Fig.  1 (left)  shows  a surface  approximated 
by  quadrilaterals. 

We  use  a quadtree  data  structure  to  keep  track  of  the  surface  subdivision 
process.  A quadrilateral  is  divided  if  it  does  not  satisfy  the  flatness  test.  Its 
children  are  subject  to  the  same  test  until  at  a certain  level  they  are  flat 
enough.  Therefore,  more  subdivisions  are  needed  at  places  where  surface 
curvature  is  high. 


§3.  Trim  Curve  Insertion 

We  assume  the  trimming  curves  are  given  in  the  parametric  space  of  the 
surfaces  which  they  delineate.  They  are  first  discretised  into  line  segments 
using  the  same  tolerance  for  the  surface.  Then  the  trimming  curve  segments 
are  inserted  into  the  quadtree  cells  by  walking  through  the  quadrilaterals  using 
adjacency  information.  The  right-hand  figure  in  Fig.  1 shows  an  example  of 
the  insertion. 

The  insertion  can  be  done  completely  in  the  parametric  domain  in  which 
the  quadrilaterals  correspond  to  rectangles  in  two  dimensions.  Starting  with  a 
vertex  of  the  trimming  segments,  we  first  find  the  rectangle  in  which  this  vertex 
is  contained.  This  can  be  done  efficiently  by  traversing  down  the  quadtree. 
By  following  the  trimming  segments,  we  can  find  the  segment  that  crosses  one 
of  the  edges  of  the  rectangle.  We  insert  a new  vertex  on  the  intersection  point 
and  then  start  the  insertion  in  the  new  rectangle. 

For  efficient  insertion  of  the  trimming  segments,  we  make  use  of  a data 
structure  that  can  quickly  find  the  neighboring  rectangle  from  the  edge  of  a 
rectangle.  We  store  in  each  rectangle  (quadrilateral)  vertex  the  pointers  of 
the  rectangles  that  use  the  vertex.  Given  an  edge  e = (p,  q),  we  collect  all 
the  rectangles  that  contain  both  p and  q,  Qe  = Qpr\Qq,  where  Qp  and  Qq 
are  the  sets  of  rectangles  associated  with  vertices  p and  q respectively.  This 
is  a local  operation.  The  number  of  elements  in  Qe  should  either  be  1 or  2. 
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Fig.  2.  Trimming  curve  insertion. 

In  the  case  of  two  elements,  one  of  them  is  the  neighboring  rectangle  we  look 
for.  When  there  is  only  one  element  in  Qe,  our  rectangle  does  not  have  a 
compatible  neighbor.  However,  if  we  go  up  the  quadtree,  at  some  level  there 
must  be  a rectangle  that  is  the  compatible  neighbor.  Then  coming  down  the 
tree,  we  find  the  leaf  rectangle  that  is  partially  neighboring  our  initial  quad. 
This  is  the  rectangle  that  the  trimming  segment  enters. 

The  time  complexity  of  the  insertion  process  is  linear  in  the  number  of 
quadrilaterals. 

Fig.  2 shows  the  insertion  of  trimming  curves  into  the  rectangles  in  para- 
metric domain.  Fig.  3 shows  the  two  cases  that  a trimming  segment  enters  a 
new  rectangle:  1)  entering  from  an  edge;  2)  entering  from  a vertex. 


Fig.  3.  Two  entering  cases. 


§4.  Initial  Triangulation 

After  the  insertion  of  the  trimming  segments,  we  have  two  kinds  of  rectangles: 
those  that  are  cut  by  the  trimming  segments  and  those  that  do  not  intersect 
with  any  part  of  the  trimming  segments.  For  each  rectangle  being  cut,  we  sort 
the  vertices  of  the  trimming  segments  inside  the  rectangle  in  counter-clockwise 
order  to  form  boundaries  of  polygons.  In  general,  there  can  be  multiple  poly- 
gons and  each  polygon  can  have  multiple  boundary  loops.  Those  cells  that  lie 
inside  the  boundary  are  triangulated  in  parametric  space.  There  is  no  short- 
age of  triangulation  algorithms  for  2-dimensional  polygonal  domains.  Here 
we  adopt  the  algorithm  from  [13].  Note  that  most  rectangles  lie  completely 
in  the  interior  of  the  trimming  boundary;  their  triangulation  can  simply  be 
done  by  triangulating  a rectangular  domain  [1],  If  a rectangle  is  cut-free  and 
lies  outside  the  trimming  boundary,  we  can  simply  ignore  it.  This  case  can  be 
decided  easily  by  testing  if  one  of  the  vertices  of  the  rectangle  is  in  the  interior 
of  the  polygon  formed  by  the  trimming  segments.  For  robustness  reasons,  we 
choose  the  centroid  of  the  rectangle  for  doing  the  test. 
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Fig.  4.  Initial  triangulation  in  parametric  space  and  in  3-space. 


An  initial  triangulation  in  3-space  is  obtained  by  evaluating  the  para- 
metric triangulation.  Fig.  4 gives  an  example  of  a surface  triangulated  in 
the  parametric  domain  (left)  and  the  corresponding  triangulation  in  3-space 
(right). 


§5.  Triangle  Shape  Improvement 

As  we  mentioned  in  Section  1,  there  are  good  reasons  to  make  triangles  that 
are  well-shaped.  In  practice,  it  is  undesirable  to  have  triangles  that  are  flat 
or  pointed.  These  are  the  triangles  that  have  one  small  angle  or  one  large 
angle.  It  is  well  known  that  Delaunay  triangulation  for  a set  of  points  in 
two  dimensions  is  optimal  in  the  sense  that  it  maximizes  the  minimum  angle 
[2],  Delaunay  triangulation  that  respect  a set  of  boundary  edges  can  be  con- 
structed. This  kind  of  triangulation  is  called  constrained  Delaunay  triangulation 
(CDT)  [12].  Chew  [3]  extended  the  definition  of  CDT  to  the  curved  surfaces 
by  replacing  the  empty  circumcircle  condition  with  the  empty  minimum  cir- 
cumsphere  condition.  Following  Chew’s  approach,  we  improve  the  shape  of 
the  triangles  by  edge  flipping  and  inserting  new  nodes  at  the  circumcenters  of 
the  ill-shaped  triangles. 

Given  a pair  of  triangles,  if  they  form  a convex  quadrilateral,  there  are  two 
choices  of  the  diagonals,  one  is  better  than  the  other  in  terms  of  the  shapes 
of  the  triangles.  By  examining  each  pair  of  adjacent  triangles  and  flipping 
their  diagonals  if  necessary,  we  can  improve  the  triangulation  locally  (see  top 
figures  of  Fig.  5).  Chew  [3]  shows  that  the  flipping  process  halts  and  it  leads 
to  constrained  Delaunay  triangulation  on  a surface. 

A CDT  is  the  best  possible  triangulation  without  introducing  new  nodes. 
To  further  improve  the  triangulation,  we  have  to  insert  new  points.  Each 
time  we  insert  a new  point,  we  do  edge  flipping  again  to  maintain  Delaunay 
triangulation.  New  points  are  inserted  at  the  circumcenters  of  the  triangles 
that  violate  the  shape  criteria.  The  reason  of  this  is  we  could  improve  the 
shape  of  several  triangles  by  introducing  one  point.  Fig.  5 illustrates  the  two 
basic  operations  used  repeatedly  for  improving  the  shape  of  the  triangles. 
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Fig.  5.  Edge  flipping  and  node  insertion. 


Fig.  6.  Example  1. 

As  we  flip  edges,  we  want  to  preserve  the  error  bound  for  the  new  trian- 
gulation. Given  a pair  of  triangles  that  are  flippable,  we  check  the  minimum 
distance  between  the  current  diagonal  and  the  new  diagonal.  The  distance 
should  be  smaller  than  the  specified  approximation  tolerance  e.  This  does 
not  guarantee  that  the  resulting  triangulation  still  satisfies  the  approximation 
tolerance.  But  since  we  are  not  moving  any  nodes  on  the  surface  and  we  in- 
sert additional  nodes  into  the  triangulation,  there  are  good  reasons  to  assume 
that  most  triangles  will  satisfy  the  tolerance.  Finally,  as  a last  step,  we  loop 
through  all  triangles  and  check  their  chord  heights.  For  those  few  triangles 
that  violate  chord  height  tolerance,  we  subdivide  them  by  adding  points  on 
their  edges.  The  checking  is  generally  expensive,  but  we  only  do  this  once 


Triangulating  Trimmed  NURBS  Surfaces 


387 


Fig.  7.  Example  2. 

for  each  triangle,  and  the  process  can  be  speeded  up  by  making  use  of  the 
quadtree  data  structure. 

Figs.  6 and  7 give  two  results  of  the  algorithm  before  and  after  shape 
improvements. 

§6.  Concluding  Remarks 

The  main  results  of  this  work  are: 

1)  a surface  triangulation  algorithm  that  guarantees  correct  mesh  topology, 

2)  an  efficient  trimming  curve  insertion  scheme, 

3)  triangle  shape  improvement  by  Delaunay  triangulation. 

We  have  only  discussed  the  problem  of  triangulating  a single  surface.  How- 
ever, in  real  world  problems,  a model  usually  consists  of  many  NURBS  surface 
patches.  The  triangulation  between  two  neighboring  surfaces  have  to  be  com- 
patible. There  should  be  a post-processing  step  that  stitches  the  triangulation 
of  different  surfaces.  This  can  be  done  with  a kd-tree  data  structure,  which 
facilities  locating  nearest  nodes  in  3-space  quickly.  Therefore,  we  can  propa- 
gate a node  on  the  boundary  of  one  surface  to  the  boundary  of  its  neighboring 
surface. 
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